function [Gamma,Su,resid,Par] = BvarGibbs_parameters(X,p,tau);

%% set tau = 0 for uninformative prior
% 
% if nargin<3
%     nScarto = 0;
%     nGibbs = 300;
% end;

[T,N] = size(X);
iRW = ones(1,N);
% Compute the parameters for the dummies
for jn = 1:N
    
    z = [ones(T-p,1)];
    
    for jp = 1:p
        z = [z X(p-jp+1:end-jp,jn)];
        y = X(p+1:end,jn);
        SS0(jn) = sqrt((y'*y-y'*z*inv(z'*z)*z'*y)/T);
    end;
    
end;

% Create matrix of regressors and dependent variable
Z = [];

for jp = 1:p
    Z = [Z X(p+1-jp:end-jp,:)];
end

Z = [Z ones(size(Z(:,1)))];

Y = X(p+1:end,:);


% Construct dummy for Litterman prior
Yrw = tau*[diag(SS0.*iRW);zeros(N*(p-1),N)];
Zrw = tau*[diag(kron(1:p,SS0)) zeros(N*p,1)];

% Construct dummy for the constant
Ycs = 1e-7*[zeros(1,N)];
Zcs = 1e-7*[zeros(1,N*p) 1];


% Construct dummies for prior on covariance matrix of residual;
Ycv = diag(SS0);
Zcv = zeros(N,N*p+1);

% put together all the information
Ypr = [Yrw; Ycv; Ycs]; Zpr = [Zrw; Zcv; Zcs];

Tpr = size(Ypr,1);

% Compute the posterior
ZZinv = inv(Zpr'*Zpr + Z'*Z);
ZY    = Zpr'*Ypr + Z'*Y;

Gamma = ZZinv*ZY;             % inv(X'X)*(X'Y) but with partitioned matrices

e  = [Y; Ypr]-[Z; Zpr]*Gamma;
resid=[Y]-[Z]*Gamma;
Su = 1/(T+Tpr)*e'*e;

k = size(Gamma,1);

CSuInv = chol(inv(Su*(T+Tpr)));
CZZinv = cholred(ZZinv);

jg = 1;


Par.CSuInv=CSuInv;
Par.CZZinv=CZZinv;
Par.Gamma=Gamma;
Par.T=T;
Par.Tpr=Tpr;
Par.k=k;
Par.N=N;
Par.p=p;

%% THIS IS THE LOOP FOR THE ESTIMATION OF THE PARAMETERS
% while jg <= nGibbs
%     
%     U    = randn(T+Tpr-k+2,N)*CSuInv;                % Inverted Wishart and Variance
%     Su_temp  = inv(U'*U);
%         
%     temp  = randn(size(Gamma));                       % Compute the betas
%     Gamma_temp = Gamma + CZZinv'*temp*chol(Su_temp);  % Instead of using KRONECHER PRODUCT!!!!!!
%     
%     AA=zeros(N*(p+1));
%     AA(1:N,1:N*p)=Gamma_temp(1:end-1,:)';
%     AA(N+1:N*p,1:N*(p-1)) = eye(N*(p-1));
%     AA(end-N+1:end,end-N+1:end)=eye(N);
%     AA(1:N,end-N+1:end)=eye(N);
%     
%     
%     Gamma_g(:,:,jg) = Gamma_temp;
%     Su_g(:,:,jg)=Su_temp;
%     jg = jg+1;
%     
% end;
% 
% Su_g = Su_g(:,:,nScarto+1:end);
% Gamma_g = Gamma_g(:,:,nScarto+1:end);


function C = cholred(S); 
% Computes Choleski factor using EIGENVALUES

[v,d] = eig((S+S')/2);

d = diag(real(d));

warning off

J = (d>1e-12);
C = zeros(size(S));

C(J,:) = (v(:,J)*(diag(d(J)))^(1/2))';

